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We use the canonical Hartree-Fock-Bogoliubov basis to implement a completely self-consistent 
| quasiparticle-random-phase approximation with arbitrary Skyrme energy density functionals and 

density-dependent pairing functionals. The point of the approach is to accurately describe mul- 
tipole strength functions in spherical even-even nuclei, including weakly-bound drip-line systems. 
We describe the method and carefully test its accuracy, particularly in handling spurious modes. 
To illustrate our approach, we calculate isoscalar and isovector monopole, dipole, and quadrupole 
strength functions in several Sn isotopes, both in the stable region and at the drip lines. 
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I. INTRODUCTION 

The study of nuclei far from stability is an increasingly important part of nuclear physics 0, @ • As radioactive 
■ beams allow more experiments on these nuclei, theoretical modeling is changing in significant ways. New ideas and 
progress in computer technology have allowed nuclear theorists to understand bits and pieces of nuclear structure 
. quantitatively |3| . Short-lived exotic nuclei offer unique tests of those aspects of our developing many-body theories 
that depend on neutron excess The major challenge is to predict or describe in detail exotic new properties of 
. nuclei far from the stability valley, and to understand the origins of these properties. 

For medium-mass and heavy nuclei, an important goal is obtaining a universal energy-density functional, which 
will be able to describe static and dynamic properties of finite nuclei and extended nucleonic matter with arbitrary 
neutron-to-proton ratio. Self-consistent methods based on density-functional theory are alreadysophisticated enough 
to allow precise analysis of ground-state properties (e.g. binding energies) in heavy nuclei 0,0,13- They can also help 
describe nuclear decays and excited states. Their predictions for collective excitations as we approach the neutron 
drip line are especially interesting. What happens to low- and high-frequency multipole modes when the neutron 
excess is unusually large? 

To address these questions we use the quasiparticle random- phase approximation (QRPA), a powerful tool for 
understanding both low-lying vibrational states and giant resonances jg. The QRPA is a microscopic approach that 
is nevertheless simple enough to allow "no-core" calculations. The approximation, which should be good for collective 
vibrations as long as their amplitudes are small, is especially effective in conjunction with Skyrme energy functionals. 
Our work is a part of a broad program to test and improve these functionals, which thus far have been fitted mainly to 
ground-state observables, by applying them to collective excitations, particularly near the drip line. This paper lays 
out our approach and evaluates its accuracy. For these purposes we restrict ourselves to a single Skyrme functional, 
SkM*. A forthcoming study will examine the performance of Skyrme functionals more generally. 

The QRPA is a standard method for describing collective excitations in open-shell superconducting nuclei with 
stable mean-field solutions, either spherical or deformed. What is not standard, and at the same time is extremely 
important for weakly bound nuclei, is the treatment of the particle continuum. Continuum extensions of the random 
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phase approximation (RPA) or QRPA are usually carried out in coordinate space, facilitating treatment of decay 
channels and guaranteeing correct asymptotics. Surprisingly, as we discuss below, the rich literature on the RPA and 
QRPA, which includes many coordinate-space calculations, contains few treatments of the continuum that exploit the 
entire Skyrme functional in a fully self-consistent way. 

To avoid confusion, we state what we mean by a fully self-consistent RPA or QRPA calculation. First, the underlying 
mean-field calculation must be self-consistent in the usual sense. Next, the residual interaction used in the RPA or 
QRPA must be derived from the same force or energy functional that determines the mean field. An important 
consequence of this condition, and of other more detailed technical conditions discussed below, is that spurious 
excitations arising from symmetry breaking by the mean field have zero or nearly zero energy, leaving the physical 
intrinsic excitations completely uncontaminated by spurious motion. Finally, energy- weighted sum rules must be 
satisfied to high accuracy. We elaborate on these requirements below; Refs. [3, Hj] discuss ways in which RPA 
calculations commonly violate them. 

The literature applying RPA or QRPA to nuclear structure is huge, and a complete review is beyond the scope 
of our paper. We do, however, present an overview of the studies that are related in one way or another to nuclear 
density functionals, self consistency, pairing, and the key issue of the particle continuum. 

The standard version of QRPA, the so-called matrix formulation, is carried out in the configuration space [T2I 
of single-quasiparticle states. A number of papers treat collective states in spherical nuclei in the Skyrme-RPA and 
QRPA matrix formulation (see Refs. (RJUjJ and references cited therein), in which the positive-energy continuum is 
discretized, e.g. by solving the Hartree-Fock-Bogoliubov (HFB) and QRPA equations in a harmonic-oscillator single- 
particle basis. Within this group, the first fully self-consistent calculations that properly account for continuum effects 
are those of Refs. 0], in which the localized canonical basis of coordinate-space HFB is used to calculate beta- 
decay rates of neutron-rich r-process nuclei and Gamow- Teller strength distributions. Recently, fully self-consistent 
HFB+QRPA calculations have also been carried out with the finite-range Gogny force 01 • Unlike many previous 
Gogny+HFB studies that employed a harmonic oscillator basis, Ref. solves the HFB equations in the eigenbasis 
of a Woods-Saxon potential, the particle continuum of which is discretized by enclosing the system in a box. 

Coordinate-space Green's functions as a method of implementing the RPA through linear response were first 
used in Ref. |18| and subsequently applied to the description of low- and high-energy nuclear modes (see, e.g., 
Refs. fl^ l2al^ 1^ 1^ 1^ 1^ 1^ l28| V Many of those calculations are not realistic enough, however, because 
they ignore the spin-orbit and Coulomb residual interactions in the RPA [ToL ITl| . Coordinate-space Green's-function 
QRPA was studied in Ref. [2!j, in the BCS approximation, with a phenomenological Woods-Saxon average potential. 
Coordinate-space HFB+QRPA for spherical nuclei was formulated in Refs. l3p . Enll32L l33l] and applied to excitations 
of neutron-rich nuclei. As in (2^, the Hartree-Fock (HF) field in Refs. |30ll3l| was approximated by a Woods-Saxon 
potential. While the calculations of Refs. j32, are based on Skyrme-HFB fields, they violate full self consistency 
by replacing the residual velocity-dependent terms of the Skyrme force by the Landau-Migdal force in the QRPA, 
and neglecting spin-spin, spin-orbit, and Coulomb residual interactions entirely. Within this approach, extensive 
Skyrme-HF+BCS QRPA calculations of El-strength in neutron-rich nuclei were carried out in Refs. |34Ll35| . 

An alternative coordinate-representation approach, also based on Green's functions, was formulated in Refs. [3rl 
l37| within Migdal's finite-Fermi-systems theory. Most of practical applications of this method, however, involve 
approximations that break self consistency in one way or another, including the use of highly truncated pairing 
spaces, different interactions in HFB and QRPA, and the so-called diagonal pairing approximation [39 . l38l l39l l4Cl l4ll 
l42l l43l liiL l4p . P roperties of excited states and strength functions have also been investigated within the relativistic 
RPAginillll m, HI HI or QRPA 0, H- The QRPA work employs the matrix formulation and is fully 
self-consistent, since it uses the same Lagrangian in the relativistic Hartree-Bogoliubov calculation of the ground state 
and in the QRPA matrix equations, which are solved in the canonical basis. 

At the present, no fully self-consistent continuum HFB+QRPA calculations exist in deformed nuclei. Refs. |55ll56| 
studied giant resonances in deformed nuclei within time-dependent HF theory, formulated in coordinate space with 
a complex absorbing boundary condition imposed . Sy mmetry-unrestricted RPA calculations, with no pairing, were 
carried out in Ref. [53 in a "mixed representation" |58j | on a Cartesian mesh in a box, while Ref. 59] contains examples 
of BCS+QRPA calculations in the single-particle basis of a deformed Woods-Saxon potential. 

The work described in this paper is fully self-consistent: among other things we use precisely the same interaction 
in the HFB and QRPA calculations so as to preserve the small- amplitude limit of time-dependent HFB. We formulate 
the QRPA in the canonical eigenbasis of the one-body particle-density matrix 60] which is calculated in the coordinate 
representation in a large spherical box. As mentioned above, the canonical basis has been used previously to study (3 
decay and Gamow- Teller strength [l5l Il6j|: its use in charge-conserving modes near the drip line is more challenging, 
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however, because of the existence of spurious states in the monopole and dipole channels. 1 These zero-energy modes 
can mix with physical states unless the QRPA equations are solved with high accuracy. A less precise implementation 
of our approach was used to calculate neutrino-nucleus cross sections in 208 Pb in Ref. [6^ |. 

This paper is organized as follows. SectionlTTlbelow presents our approach. In Sec. lIIII we check the QRPA solutions 
carefully, focusing on spurious modes. Section llVI contains the main conclusions of our work. Mathematical details 
are in two appendices, the first of which is on the QRPA equations and the second on calculating the derivatives of 
the Skyrme functionals that enter the formalism. 



II. METHOD 



Our first step in the self-consistent treatment of excitations is to solve the spherical HFB equations in coordinate 
space (without mixin g ne utron and proton quasiparticle wave functions |63|). with the method developed in Ref. |64j 
(see also Refs. [60H65ll6^ |'). We can use arbitrary Skyrme functionals in the particle-hole and pairing (particle-particle) 
channels. 

We modify the code used in Refs. [60l 163. l65| so that it solves the HFB equations with higher accuracy, which 
we need because the QRPA uses all the single-quasiparticle states produced by the HFB equations, even those that 
are essentially unoccupied. Our modifications are: (i) the use of quadruple precision (though in solving the QRPA 
equations we use double precision); (ii) a smaller discretization length (0.05 fm); and (iii) a high quasiparticle-energy 
cutoff (200 MeV) and a maximum angular momentum j max =15/2 (N < 82) or 21/2 (N > 82). In a 20 fm box, 
this cutoff corresponds to 200-300 quasiparticle states for each kind of nucleon. We include all these quasiparticle 
states in the HFB calculation because a very large energy cutoff is essential for the accuracy of self-consistent QRPA 
calculations ^(j. Hence, the effective pairing window in our HFB calculations is also very large, with the pairing 
functional fitted to experimental pairing gaps extracted as in Ref. 67] from the measured odd-even mass differences 
in several Sn, Ni, and Ca isotopes. 

Next, we construct the canonical basis, the eigenstates of single-particle density matrix p. To avoid poor accuracy 
(see Ref. js^l ) m the wave functions of the nearly empty canonical particle states, we do not diagonalize p directly in 
coordinate space. Instead we construct an intermediate basis by orthonormalizing a set of functions {ip^(r) + tp%(r)\, 
where <fi(r) and <P2 ( r ) are the upper and lower components of the quasiparticle wave function with energy [64|. 
We use the density matrix in coordinate space to calculate the matrix in this basis, which we then diagonalize to 
obtain the canonical states. The reason for using the sum of fi(r) and ^ ( r ) is that solutions of the HFB equations 
expressed in the canonical basis (Eqs. (4.14) of Ref. |60j| ') are, in the new basis, guaranteed to be numerically consistent 
with those of the original HFB problem. This is because the configuration space is the same in both cases, independent 
of the pairing cutoff (see Ref. (H3 for a discussion relevant to this point). Without pairing, when either tp^ (r) or 
iP2 (r) is equal to zero, our method is equivalent to taking a certain number of HF states, including many unoccupied 
states. 

In the canonical basis, the HFB+QRPA equations have a form almost identical to that of the BCS+QRPA ap- 
proximation, the only difference being the presence of off-diagonal terms in the single-quasiparticle energies. The 
QRPA+HFB formalism employs more pairing matrix elements than the QRPA+BCS, however. 

As noted already, full self consistency requires the use of the same interaction in the QRPA as in the HFB approx- 
imation. More specifically, this means that the matrix elements that enter the QRPA equation are related to second 
derivatives of a mean-field energy functional. We describe the densities and the form of the functional carefully in the 
appendices. But we must meet other conditions as well for QRPA calculation to be self-consistent. Essentially all the 
single-particle or quasiparticle states produced by the HFB calculation must be used in the space of two-quasiparticlc 
QRPA excitations. This requirement is rather stringent, so we truncate the two-quasiparticle space at several levels 
and check for convergence of the QRPA solution. First we omit canonical-basis wave functions that have occupation 
probabilities v\ less than some small v^ rit , (or HF energies greater than some e cr it if there is no pairing). Then we 
exclude from the QRPA pairs of canonical states for which the occupation probabilities are both larger than 1 — v^ lit . 
This second cut is based on the assumption that two-particle transfer modes are not strongly coupled to particle-hole 
excitations. In addition, if the factors containing Ui and Vi in the QRPA equation — see Eqs. (| Al 2|) and l|A13Jl — are 
very small, in practice smaller than 0.01, then we set the corresponding matrix elements equal to zero. This does not 
affect the size of the QRPA space, but significantly speeds up the calculations. For good performance we diagonalize 
QRPA-Hamiltonian matrices of order 20,000 x 20,000 in neutron-rich Sn isotopes. 



As far a s w e know, the only application of the canonical basis to charge-conserving modes near the drip line is in the relativistic QRPA, 
see e.g. ISSIrJU . 
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Having solved the QRPA equations, we can then calculate the strength function 

for the multipole operator Fjm ■ The smoothing width 7 is supposed to be large enough to remove spurious oscillations 
in Sj(E) associated with a finite box radius i?box |55ll69j - A reasonable form, based on a single-particle estimate, for 
the smoothing width (App. B of Ref. [(59)), is 



~f(E) = { R^V ri2( 2m An) ' E >-^n (2 ) 

[ 0.1 MeV, E < -A n 

where A n is the neutron Fermi level and m is the nucleon mass. In deriving Eq. (J2J) we assumed that the single-proton 
continuum is effectively shifted up several MeV by the Coulomb barrier. In other words, we associate the threshold 
energy with the neutron Fermi level. 

In all the tests below, we use the Skyrme functional SkM* Jjjj and a volume pairing functional [TJ (C^(poo) a 
constant in Eq. l)B9jO . The pairing parameter in Eq. (|B19(I is Vq = —77.5 MeV fm 3 . Usually we work in a box of 
radius 20 fm, though we vary this radius below to see its effects. In several tests we examine the weakly bound nucleus 
174 Sn, which is very close to the two-neutron drip line. In this system, the protons are unpaired and the neutrons 
paired (with A n =1.016MeV) in the HFB ground state. 



III. ACCURACY OF SOLUTIONS 

Benchmark tests of the HFB part of our calculations are reported in Ref. jl3 ■ Since the accuracy of the canonical 
wave functions, in which the QRPA calculations are carried out, strongly affects the quality of results (in particular 
QRPA self consistency), we take special care to compute them precisely. As discussed in Sec.|n] we obtain canonical 
states by diagonalizing the single-particle density matrix p represented in the orthonormalized set of functions {ip^(r) + 
if>2 (?")}• The accuracy of this method is illustrated in Fig. 2] which plots the quasiparticle energies £^ hock , obtained by 
diagonalizing the HFB Hamiltonian in the canonical basis (Eq. (4.20) of Ref. |6(|), versus the quasiparticle energies 
E^ obtained by solving the HFB differential equations directly in coordinate space (Eq. (4.10) of Ref. [(5(j). Two 
sets of canonical states are used: (i) those obtained through the procedure outlined above (dotted line) and (ii) those 
obtained in the standard way by diagonalizing the density matrix p(r, r') in discretized coordinate space (Eq. (3.24a) 
of Ref. solid line). If the canonical basis is precisely determined, E'^ loc]t =E fl and the two sets of £;^ heck coincide. 
Within the standard approach, however, the high canonical energies deviate visibly from their HFB counterparts, i.e., 
the accuracy of the underlying canonical wave functions is poor. On the other hand, the quasiparticle energies and 
canonical wave functions calculated within the modified approach introduced above are as accurate as the original 
solutions to the HFB equations, even for high-lying nearly-empty states. (See also Sec. VI. D of Ref. [73j for a discussion 
relevant to this point.) 

Having examined the canonical basis, we turn to the accuracy of the QRPA part of the calculation. To test it, we 
first consider solutions related to symmetries. If a Hamiltonian is invariant under a symmetry operator P and the 
HFB state \*&) spontaneously breaks the symmetry, then e iaP \9), with a an arbitrary c- number, is degenerate with 
the state |^). The QRPA equations have a spurious solution at zero energy associated with the symmetry breaking 
while all other solutions are free of the spurious motion. This property is important for strength functions, 
and gives us a way of testing the calculations. Since our QRPA equations, which assume spherical symmetry, are 
based on mean fields that include pairing and are localized in space, there appear spurious states associated with 
particle-number nonconservation (proton and/or neutron; + channel) and center-of-mass motion (1~ channel). These 
two cases are discussed below in Sec. 1111 Al and Sec. 1111 Bl 



A. The isoscalar mode 



In addition to the spurious state associated with nonconservation of particle- number by the HFB, the + channel 
contains the important "breathing mode". In Table [I] we display results from a run with w^ rit =10~ 12 for neutrons 
and e cr it = 150 MeV for protons, resulting in the inclusion of 310 proton quasiparticle states and the same number 
of neutron states, with angular momentum up to j = 21/2. The Table shows the QRPA energies and transition 
matrix elements of the particle-number operator. The spurious state is below 200 keV, well separated from the other 



FIG. 1: Neutron quasiparticle energies E^ oc for s 1 / 2 states in Sn, calculated by diagonalizing the HFB Hamiltonian in the 
canonical basis, versus the quasiparticle energies obtained by directly solving the HFB equations in coordinate space. Stan- 
dard (solid line, diamonds) and improved (dotted line, dots) methods are used to obtain the canonical states. See text for details. 



states, all of which have negligible "number-strength" . The nonzero number strength in the spurious state, like the 
nonzero energy of that state, is a measure of numerical error. If the space of two-quasiparticle states is smaller, with 



TABLE I: The lowest-energy excited + states in 174 Sn. The second column shows the excitation energies and the third 
column the squared matrix elements of the particle- number operator between the fcth excited state and the ground state (k=0). 



k E k (MeV) |{fc|jV|0)| 2 

1 0.171 0.120 

2 2.833 0.533 x 10~ 5 

3 3.090 0.877 x 10~ 7 

4 3.810 0.252 x 10 _s 

5 3.878 0.480 x 10~ 5 
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FIG. 2: Isoscalar + strength function in 174 Sn for (i) the single-proton energy cutoff e cr it = 100 MeV and the neutron- 
quasiparticle occupation cutoff u 2 rit = 10~ 8 (thin solid line); (ii) e cr it = 150 MeV and u 2 rit = 10~ 12 (dotted line); and 
(iii) £ C rit=200 MeV and f 2 rit = 1CP 16 (thick solid line). Results corresponding to (ii) and (hi) practically coincide. 



£crit = 100 MeV and t> 2 rit = 10 8 , the energy of the spurious state and the number strength barely change. 
Figure |2 shows the strength function Sj(E) for the isoscalar + transition operator, cf. [75| . 

£>o = ^X>'- ( 3 ) 

i=l 

We have plotted three curves with successively more quasiparticle levels (from 246 proton levels and 203 neutron 
levels to 341 proton levels and 374 neutron levels), with cutoff parameters given in the figure caption. The major 
structures in the strength function are stable. The error remaining after to the gentlest truncation is extremely small. 

The dependence of the strength function on the box size and quasiparticle cutoff is shown in Fig. [3J The upper 
part of the Figure (panels a-c) corresponds to a constant smoothing width of 7=0.5 MeV. This relatively small value 
is not sufficient to eliminate the finite-box effects but it allows us to assess the stability of the QRPA solutions as a 
function of i?box- The large structure corresponding to the giant monopole resonance (GMR) is independent of box 
size no matter what the cutoff, but increasing the number of configurations magnifies the dependence on box size of 
local fluctuations in Sj(E). The lower part of the Figure (panels d-f) are smoothed more realistically, as in Eq. 
It is gratifying to see that the resulting strength functions are practically identical, i.e., the remaining dependence on 
i?box and the cutoff is very weak. 

The energy-weighted sum rule (EWSR) for the isoscalar + mode [75|| is given by 

^^|(fc|^ 00 |0)| 2 = 2----^(r 2 ), (4) 

k 

where the expectation value is evaluated in the HFB ground state. This sum rule provides a stringent test of self 
consistency in the QRPA. In 174 Sn, the right-hand side of Eq. @ is 35215 e 2 MeV fm 4 and the left-hand side 34985±15 
e 2 MeV fm 4 for all of the calculations of Fig. |3 the QRPA strength essentially exhausts the sum rule. (The QRPA 
values of the EWSR in this paper are obtained by summing up to Ek — 50 MeV. ) 



B. The isoscalar 1 mode 



The 1 channel, home of the giant dipole resonance, the isoscalar squeezing resonance, and as yet incompletely 
understood low-energy peaks in neutron-rich nuclei (sometimes associated with skin excitations), has a spurious 
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E(MeV) 

FIG. 3: Isoscalar + strength function in 174 Sn for the box radii: i?box=20fm (solid line) and J?box=25fm 
(dotted line). In (a), (b), and (c) the smoothing- width parameter 7 is 0.5 MeV for all energies, while 
in (d), (e), and (f) y(E) is given by Eq. (|5J. We use the same three sets of cutoff conditions as 
in Fig. [5] namely (i) in parts (a) and (d), (ii) in parts (b) and (e), and (iii) in parts (c) and (f). 



isoscalar mode associated with center-of-mass motion that can seriously compromise the low-energy spectrum if not 
handled with extreme care. We test the ability of our QRPA to do so in 100 Sn, 120 Sn, 174 Sn, and 176 Sn. (The nuclei 
100 Sn and 176 Sn are the two-proton and two- neutron drip- line systems predicted by the HFB calculation with SkM*. 




FIG. 4: Isoscalar 1 strength function in 100,120,174,176 gn 

for the corrected dipole operator in Eq. © (solid line) 
and the uncorrected operator in Eq. J5J (dotted line). The cutoff e C rit is 140 MeV and v* lit is 3 x 10" 12 . The self 
consistency of our calculations makes the solid and dotted curves coincide nearly exactly over the whole energy range. 
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Neither nucleus has any static pairing, i.e., A n =A p =0.) In the following calculations, we take e C rit = 140 MeV for 
the protons and v% lit = 9x 10~ 12 for the neutrons. As discussed above, smoothed strength functions are practically 
independent of small changes in the cutoff. They are also independent of the cutoff in quasiparticle angular momentum 
provided we include all states with j<15/2. 

Figure m shows the predicted isoscalar dipole strength function for 100 . 120 < 174 > 176 s r L For the transition operator, we 
use 



Fim = ^j-J2rfY 1M {tti 



A 



i=l 



(5) 



and the corrected operator 



TT'COr 

#1M 



7 A 

eZ % - 

i=l 



(6) 



to remove as completely as possible residual pieces of the spurious state from the physical states ^(| • The fact that 
the strength functions produced by these two operators — displayed in Fig.0] — coincide so closely shows the extreme 
accuracy of our QRPA solutions; they are uncontaminated by spurious motion even without the operator correction. 
The spurious-state energies -©spurious are 0.964 MeV for 100 Sn and 0.713 MeV for 120 Sn, and the energies of the first 
physical excited states are 7.958 MeV for 100 Sn and 7.729 MeV for 120 Sn. In 174 Sn ( 176 Sn), Spurious is 0.319 MeV 
(0.349 MeV) and the first physical state is at 3.485 MeV (2.710 MeV), lower than in the more stable isotopes. Pairing 
correlations do not affect accuracy; the neutrons in 120 Sn and 174 Sn are paired, while those in 100 Sn and 176 Sn are 
not. 

We display the fine structure of the isoscalar 1~ strength functions in 120 Sn and Sn in Fig.0 which also illustrates 
the dependence of the results on i?box- The dependence is consistent with that of Fig.[3]for the isoscalar + strength; 
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FIG. 5: Isoscalar 1~ strength function in 120 Sn (left) and 174 Sn (right) for two box 

radii: _Rbox=20fm (solid line) and -Rbox=25fm (dotted line). In (a) and (b) the smoothing- 

width parameter is constant (7=0.5 MeV), while in (c) and (d) 7(-E) is given by Eq. 
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the low-amplitude fluctuations in Sj(E) that are unstable as a function of i?box disappear, and the smoothed strength 
function depends only weakly on i?box- In 120 Sn, the two sharp peaks below 10 MeV correspond to discrete states 
while the broad maxima centered around 15 MeV and 27 MeV are in the continuum, well above neutron-emission 
threshold. A similar three-peaked structure emerges in 174 Sn, though most of the strength there is concentrated in the 
low-energy peak at E f» 4 MeV. Fig. 0] shows (as we will discuss in our forthcoming paper [76|) that the appearance 
of the low-energy isoscalar dipole strength is a real and dramatic feature of neutron-rich dripline nuclei [23, [Za] ■ The 
EWSR for the isoscalar 1~ mode [75| is 

EE^<^mI°>| 2 = (V) - f (rr) • (7) 

ft M V 7 

In 174 Sn, the right-hand side is 403310 e 2 MeV fm 6 , while the left-hand side is 400200 e 2 MeV fm 6 . For 176 Sn, the 
corresponding numbers are 406576 e 2 MeV fm 6 and 407100 e 2 MeV fm 6 . This level of agreement is very good. 



C. The isovector + and isoscalar, isovector 2 + modes 

Figure El displays strength functions for the + and 2 + channels in 120 Sn and 174 Sn. (We discussed the isoscalar + 
mode above to illustrate the accuracy of our solutions, but include it here as well for completeness.) The calculations 
show the appearance of low-energy + strength — both isovector and isoscalar — and low-energy isovector 2 + strength 
in 174 Sn, though in none of these instances is the phenomenon quite as dramatic as in the isoscalar l - channel. 

The EWSR for the isoscalar 2 + transition operator, 

Z A 

F2M = e-Y j r 2 l Y 2M (n i ), (8) 




FIG. 6: Isoscalar and isovector strength functions for (a) the + channel of 120 Sn, (b) the + channel of 174 Sn, 
(c) the 2+ channel of 120 Sn, and (d) the 2+ channel of 174 Sn. The cutoff e crit is 150 MeV and v* lit is 1(T 12 . 
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can be written as 

EE^i< fc i^'i°>i 2 = r— w 

fc M 

The sum rule is obeyed as well in the 2 + isoscalar channel as in the + and 1~ channels, the only difference being 
that one needs to include quasiparticle states with j > 15/2 for 174 Sn. For 120 Sn ( 174 Sn) from Fig.EJ the EWSR is 
37222 (34971) e 2 MeVfm 4 while the QRPA value is 37030 (35010) e 2 MeVfm 4 . 

While on the topic of the sum rule, we display in Table ITTI the j max -dependence of the EWSR for several channels 
in 150 Sn, with i?box = 25 fm. By taking j ma x = 19/2 we appear to obtain essentially the entire strength in all three 
cases. 



TABLE II: The j max -dependence of isoscalar EWSR for 150 Sn. i? box is 25 fm. 



r ,r 


units 


Jmax — 19/2 


jmax — 25/2 


IS 0+ 


[e 1 MeV fm 4 ) 


35731 


35633 


is i- 


(e 2 MeV fm 6 ) 
(e 2 MeV fm 4 ) 


361686 


353936 


IS 2+ 


35542 


35445 



IV. CONCLUSION 



In this work we have reported on the development and detailed testing of a fully self-consistent Skyrme-QRPA 
framework that employs the canonical HFB basis. The method can be used to calculate strength distributions in any 
spin-isospin channel and in any spherical even-even nucleus. A good calculation requires a large single-quasiparticlc 
space. Our results show that our space is large enough in nuclei as heavy as the Sn isotopes. 

We are currently investigating the predictions of a range of Skyrme functionals across the Ca, Ni, and Sn isotope 
chains. The initial results presented here point to increases in low-lying strength at the neutron drip line, particularly 
in the isoscalar-dipole channel. In a forthcoming paper j7|| we will report on the robustness of these effects, on the 
physics underlying them, on their variation with atomic mass and number, and on their implications for the future 
of Skyrme functionals. 
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APPENDIX A: QRPA EQUATION 

The QRPA equations are the small-oscillations limit of the time-dependent Hartree-Fock-Bogoliubov approximation, 
see, e.g., Hjjj. In the canonical basis the most general equations take the form 



E 

L<L' 



AkK'.LL' BkK',LL' \ ( X Lr , \ „ I X\ 

rk 



B * _ A * j i Jk V )=Ek[ yV , K<K', (Al) 

KK'.LL' KK' ,LL' J \ r LL> / \ Y KK' 
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Akk\ll' — EklSk'L' — Ek'lSkl' — Ekl'Sk'l + Ek>l'Skl 

-V^ h LR , L ,UL>V L UKVK' + V^ iRL ,UL'VLUK'VK 
+ V^ h L , R , L U L VL'UKVK' - V^, L , Rl ULVL'UK'VK 

-VlIk'k VlVl ' Vk ' Vk - 

-VlL'kk ,VlVl ' UkVk ' + VujK' K VlVl ' Uk ' Vk 
-V r P j^ Ll ,ukUk'U L >v l + V r p r , l , l u k u k >u l Vl' 
-Vll^, R u L >v L v K 'VK + V L P l RlR u L v L ,v K ,v K 

-V r p r , l , l u k v K 'Ulu L ' +V R p R h L , L u K >v K u L u L -, (A2) 
Bkk'xu = V p , L , RL u L 'V L u K 'V K - V ph 



V W'T kt.UlVL'Uk'Vk + V ph r f>, T ,U L V L 'UKVK' 



K' LKL' ^ i\ i KLK'V 

+ V PP rLL ,V L VL'UkUk> + Vl P LRRl VKVK'U L U L ' 

+ V R P R h L , L UL>V L UKUK> - V R P R h LL ,ULVL>UKUK> 

+Vl?l R , R u K ,v K u L u L , - Vl P l KRl u R v R ,u L u L , 
+V R p RLL ,v L v L ,u R ,v R - V K v R ] l LL ,v L v L ,u R v R , 

+Vl P l RRI v K v K >u Ll v L - vII^ rR ,v k vk>ulVl> , (A3) 



S 2 E[p : K, K* 
SpK'KSpL'L 



c = ^rr;::. . ( A4 ) 



_ 5 2 E[p, K ,K*} 

Vk ' kl ' l - 5k k , k 5k l ^ (A5) 



x>3plh _ 8 E[p, K, K j _ 7 ip 3 h * , 

v K'KL'L — ~T~i 77 — V LL'K'K> \ A0 ) 
dh R , R dp L 



where K and L are single-particle indices for the canonical basis, and the states are assumed to be ordered. The 
symbol K refers to the conjugate partner of K, uk and vk come from the BCS transformation associated with the 
canonical basis, and the Ekl are the one-quasiparticle matrix elements of the HFB Hamiltonian (cf. Eq. (4.14b) 
of Ref. jSfl|). A| L , and Y^ L , are the forward and backward amplitudes of the QRPA solution k, and Ek is the 
corresponding excitation energy. E[p, n, k*] is the energy functional (see Add. iBlfor an explicit definition) and p and 
k are the density matrix and pairing tensor. After taking the functional derivatives, we replace p and k by their HFB 
solutions, in complete analogy with an ordinary Taylor-series expansion. 
To write the equations in coupled form, we introduce the notation 

K = {n^l^j^m^) = (pm^), L = {vm v ), (A7) 

where (nljm) denote spherical quantum numbers. Using (i) rotational, time-reversal, and parity symmetries of the 
HFB state, (ii) the conjugate single-particle state 2 

\K) = \W^) = (-Y^\p -m^), (A8) 



2 The conjugate state in our HFB code is slightly different: \fJ,m) = (— ) J m+l \fi — m). This definition follows from Eq. (3.24b) of Ref. |60H 
and the single-particle wave function ipfimir) = Rn( r ) Xw 2t r §°"li m ) \ a }> where -R^(r) and \a = ±1/2) are real radial and spin 

wave functions. Thus, we multiply all HFB by (— ) ; in the QRPA calculations. 
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and (iii) the relations 



X* KK . = {3^ imM M k )Xl, ]Jk x { f q ^-^ (A9) 



Y, 



KK' 



{-y^i-y^^-m^-m^J^Y^ x Y' otherwise (Al0) 



fe — f—\U-m u (_\3 u ,'-m u ,i /„• _ m „• . _ ™ . I T. A/f. WI* . v J V^j H — H' , 



with Jfc the angular momentum of the state k and the factor \/2 for convenience |l2j| . one can rewrite the QRPA 
equation as 

/ ( jWu-v* JW*.^ \ ( K v , ]Jk \ =E ( V < fi', (ah) 



At, in,' 



( w; 
,/ /. 



+G(fj,fi vv ; J^iu^u^UyUv' + VvVisVpiVfi) 

-(-)^'+>- J *F(^/xW; Jk){u,MV v u v >v^ + Up'V^UvVp) 

+(-y»+iu'- J >*H((j,'(U'p'; JkjivvVy'U^v^ + u ft 'iyu„'v) 
-H*(w'fx/j,'; JkjiuiiUp'UvfVv + u„?v 

+(-y"+iv'- J *H*(i/'i/(iti']J k )(u IJ ,u tt ,u u v„> + u u .v v v ll .v^) }, (A12) 

= /-, 1 x /-, 1 x {-G{pn'vv';J k ){u ll ,u ll VuV v <+u v u v ,v ll >v ll ) 
V 1 + <W v 1 + <W 

-{-y v+jv '~ Jk F(niJ,'v'v : Jk^U/iUvVu'Vv' + u^Uv'VvV^) 

_)_(_ y»+3»'+3»+0it' J?(n' [j,]/ V ; J^Ufjt'UuVviVfj, + UfiUv'Vfi'Vu) 

+H{nn'vv'\ JkjivuVu'UpiVp + u^v^UvUv*) 
-H*(pp'fj,fjf; Jfc)(w Al v^u v iv v + UuV^UftUft') 

+{-y" +j "'~ Jk H* {v 'i/fifi' '; Jk)(vpVpiu u v v > + u u >v v UpUp>) } , (A13) 
G(p(i'vi/;Jk) = ^ (jvmvjv'mv^JkM k }(j„m„j„,m^\J k M k }V£ p K , LL , 

m^m r m v m v i 

= {W]J k \V pp \[ VV ']J k ), (A14) 
F(/i/zW;J fe ) = {j^nipi^m^\J k M k )(j v mvj v >mv'\JkMk)V^ LlR , L 

= 5>)V+^' | j V «j* | (2 J' + 1)([^'] J fc |?P h |[A] J k ), (A15) 
H{nn'w'\J k ) = ^2 {j^^j^m^\JM){j v m v j v im v '\JM)Vlf; KRl 

= V(-^+^+ 1 + , -'- J *- J ' ( j« ^ )(2j' + i)([H^|^ 3plh |[^VV fc >. 



(A16) 



We have represented the second derivatives of the energy functional E[p, k, k*] as unsymmetrized matrix elements of 
effective interactions V pp , V ph , and V 3plh . These effective interactions are given in App. IbI Although the "matrix 
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elements" are unsymmetrized, the underlying two-quasiparticle states are of course antisymmetric. As a consequence, 
A.[ l i l i'}j k ,[uu'}J h = B [^']J k \^'].J k = if Jk is odd and either fi = fi'ori/ = v' . 

The nuclear energy functional E[p, k, k*] is usually separated into particle-hole (ph) and pairing pieces (again, see 
App. B for explicit expressions). If the pairing functional, which we will call E na ,\ r \p, k, k*], depends on p then the 
derivatives of E pa i r [p, k, k*] with respect to pkk' are called pairing-rearrangement terms [13j . In the QRPA, two kinds 
of pairing-rearrangement terms can arise in general. One has particle-hole character and is included in V^ LK i L ,; the 
other affects 3-particle-l-hole (3plh) and l-particle-3-hole (lp3h) configurations and is represented by Vv!\cl,'L an< ^ 
^k'^kl'l- If the p-dependence of -Epair [p, k, k*] is linear, then the ph-type pairing-rearrangement term does not appear. 
Furthermore, the 3plh and lp3h pairing-rearrangement terms arise only for J 7 ' = + modes if the HFB state has 
J = 0. Most existing work uses a pairing functional that is linear in p, and so needs no pairing-rearrangement terms 
in J w 7^ + channels. 



APPENDIX B: INTERACTION MATRIX ELEMENTS (SECOND FUNCTIONAL DERIVATIVES) 

1. Representation of second derivatives as matrix element of effective interactions 

In this appendix, we discuss interaction matrix elements coming from E[p, k, k*], which we take to contain separate 
Skyrme (i.e. strong-force, K-independent) , Coulomb, and pairing energy functionals: 

E[p, K, K*] = BskyrmeW + #Coul[p P ] + #pair[p, K, «*], (Bl) 

where p p is the proton density matrix. The most general Skyrme energy functional in common use is given by 
E Skyimc [p] = E / M CP M Pto(r) + C? p p t o(r)Ap m (r) + C\ [p t0 (r) r tQ (r) - f m (r)] 

4=0,1*' 

+CM s 2 t0 {r) + C^s t0 (r) ■ As t0 (r) + Cj [s 00 (r) ■ T t0 (r) - J 2 m (r) 

+C? J [pt (r)V • Jto(r) + s m (r) ■ V x j t0 (r)} 1. (B2) 



(See, e.g., |63l l79l |8fj| and references therein for a general discussion.) All densities are labeled by isospin indices tt z , 
where t takes values zero and one and t z is always equal to 0. A more general theory could violate isospin at the 
single-quasiparticle level, leading to additional densities pi±i H^. We do not consider such densities here. The CI 
are the coupling constants for the effective interaction. As usual, two of them are chosen to be density dependent: 

Cf[ P oo] = AP+BPp^ (r), 

CM = At + BlpUr)- (B3) 

Here p t Q, s t o, Tto, Tto, jto, J to, and Jto are local densities and currents, which are derived from the general density 
matrices for protons and neutrons 



p 00 (ra,rcr) = p n (ra,r a ) + p p (ra,r a ), 

p w (ra,r'a') = p n (ra,r'a') - p p (ra,r'a'), (B4) 

where 

p a (ra,r'a') = ^ t()g,(r'a , )ipK(rcr)pjcK', 

KK' , neutron 

p p (ra,r'a') = ^ ^K'{r'a')i) K (r(i)pKK', (B5) 

KK' , proton 

and a — ±i labels the spin components so that, e.g., ipj((ra) is a spin component of the single-particle wave function 
associated with the state K. Defining 



PtO 



Sto(r,r') = ^2 Ptaira.r 1 a')a a : ai (B6) 



a.a'— ± 
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= (a'\cr\a) is a matrix element of the vector of Pauli spin matrices, we write the local densities and 
currents as 

Pto(r) = fho(r,r), 

s m (r) = s t0 (r,r), 

Tto(r) = V- Vpta(r,r')\ r=r >, 

T m (r) = V- V's t0 (r,r')\ r=r ,, 

3m(r) = -i(V-V')pto(r,r')\ r=r ,, 

Jto,ij(r) = -|(V - V')iS t oj(r,r')|r=r', 

J to( r ) = ^tO,iji 
ij=xyz 

Jto(r) = -|(V- V')xs tQ {r,r')\ r = r ,. (B7) 
The Coulomb energy functional is given by 

EcM = & j J J ArfV^ - \? (!) 1 / <fV #»(r), (B8) 

where we make the usual Slater approximation |8l| for the exchange term. 
For the pairing functional we take the quite general form 

E paiT [p,K, K *] = E pail [p,p,p*} = J d 3 r CP\p 0Q (r)] £ |/5 T (r)| 2 , (B9) 

r— p.n 

where the density-dependent pairing coupling constant C p [poo( r )] is an arbitrary function of poo(r). The quantity 
p T (r) is defined as 

p T (r) = —i K T (ra, r<j')o v aa , , r = proton or neutron, (BIO) 

<t<t'=± 

with 

K n (ra-,rV) = %l> K .(r' (j')%}; K (ro)K KK , , 

K K ' , neutron 

K p (rcr,rV) = i} K/ (r'a')i} K (ra)K KK ,, (Bll) 

KK' , proton 

being the standard pairing tensor in the coordinate representation. 

The second derivatives of the energy functional in Eq. (|A4|) . as the equation indicates and we've already noted, 
can be written as unsymmetrized matrix elements V^lk'L' °f an effective interaction between uncoupled pairs of 
single-particle states. The particle-hole matrix elements take the form 

V$lk>u = (KL\V^ ymc + Vf oul + V;l* h \K'L'). (B12) 

The last term contains the pairing rearrangement discussed at the end of the previous appendix. 
The effective Skyrme interaction in Eq. (|B12|) is given by 



^Slfyrmc = K + b <T • <t' + C f • f + d Q (7 ■ a' f ■ f) 5(r - r') 

+(ai + bi a ■ a' + c x f ■ f 1 + di a ■ a' f ■ t") (fc f2 S(r - r') + 6(r - r')k 2 ) 
+(a 2 + b 2 cr ■ cr' + c 2 t ■ f' + d 2 cr ■ cr' f ■ f') fc 1 * • 5(r - r') k 
+(a 3 + b 3 cr ■ cr' + c 3 r • f + d 3 cr ■ cr' r ■ 7") p% Q {r) 5(r - r') 
+ [e 3 Pw (r) (r<°> + r'(°)) + g 3 s 00 (r) ■ (cr + cr') 

+m 3 s w {r) ■ (art ' + cr'r'^)} p a m \r) S(r - r') 
+ [h Ploir) + h 3 s 2 Q0 (r) + n 3 s 2 10 (r)] p° Q 2 (r) S(r - r') 

+(a 4 + c 4 t ■ f') (cr + cr') • fet x j( r _ r ') k (B13) 
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TABLE III: Definitions of en, bi, Ci,di (i = 0, ■ ■ ■ ,3), 04, and C4 in Eq. 1HT31. 



i aj_ bi Cj di 

2A% 2lg 23| 

1 |(q-4C) I(C T -4C As ) I(CI-4Cf) l{cT-ACt) 

2 3C T +4C Ap 3C T +4C^ S 3<7f+4C 1 v 3d T + 4C As 

3 B£(a + 2)(a + l) 2Bg 2Bf 2B{ 

4 - 2iC,7 J - 2iCY J 



TABLE IV: Definitions of the coefficients appearing in the rearrangement terms in Eq. 1B13II . 



i 


&i 


fi 


9i 


hi 


m,i 


rii 


3 


2aB{ 


a(a - l)Bf 


2aB^ 


a{a - 1)Bq 


2aB{ 


a(a - l)Bi 



t = (r^ 1 ^ , ) is the vector of Pauli matrices in isospin space and 

k = — |(V — V) acting to the right, 

fe f = l(V-V') acting to the left. (B14) 

The coefficients in Eq. (IB13I) are defined in Tables ITTll and llVl Equation (|B13(1 contains the usual Skyrme- interaction 
operators, but, the energy functional l|B2|) does not necessarily correspond to a real (density-dependent) two-body 
Skyrme interaction because the matrix elements are not antisymmetrized. Compared to the case usually discussed in 
the literature, the more general functional relaxes relations that would otherwise restrict the spin-isospin structure of 
the effective interaction in Eq. (|B 13|) : see, e.g., 0] for a discussion of the increased freedom. 

The densities and currents that appear in Eq. (|B13|) come mostly from rearrangement terms and take the values 
given by the HFB ground state. The isoscalar and isovector spin densities s t o(r) vanish when the HFB ground state 
is time-reversal invariant or spherical as assumed here. The terms containing them will therefore not appear in the 
expressions for the matrix elements of the effective interaction for such states given below. 

The effective Coulomb interaction in Eq. I|B12|I . acting between protons, is given by 



^ = r-^-jr ~ t (z) ' P^) 8(r - r'). (BIS) 



TV 

Finally, the ph-type pairing-rearrangement terms in Eq. (|B 12|) come from an effective interaction 

rC s P h _ d 2 CP[p o(r)} 

t— p,n 



jgph = a wpoop E |Mr)|2 §(r , r)< (B16) 



The second derivatives with respect to k, k* also can be written as unsymmetrized matrix elements of effective 
interactions, this time in the particle-particle channel. The particle-particle effective interaction entering the matrix 
elements 

n%LL. = (KK'\V;«r\LL') (B17) 
is obtained from Eq. (|B9|I through Eq. (|A5|) as 

V^ pp = C*boo(r)] (3 - a ■ a' - f ■ t> - a ■ a' f ■ f') S(r - r'). (B18) 
In the numerical calculations of this paper, we use a volume pairing-energy functional, i.e., 

= -Vq = const. (B19) 

Last of all are the mixed functional derivatives involving both p and k (or p) in Eq. (|A6|I . They also can be written 
as the unsymmetrized matrix elements of an effective interaction: 

y$KVL = (L'K'\V;l 3plh \LT(K)), (B20) 
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where T(K) denotes the time-reversed state of K, and the 3plh effective interaction itself is 

Vpl 3plh = ^[r) 01 [ ^ (r) (1 " <] + Mr) (1 + T ^ )] ^ " r) ' (B21) 

where t' z acts on the single-particle states K' and T(K) in Eq. I|B20|) . and the eigenvalues 1 and —1 are assigned to 
the neutron and proton, respectively. 



2. Calculation of matrix elements 

To calculate the coupled matrix elements in Eqs. I|A14|) - (|A16|) . we use an intermediate LS scheme: 

(\fin']J k \V\[vt/]J k ) 




Y, hhhh'LL'S 2 I 1/2 1/2 5 

LL'S { jit V J k 

xdy^LS; J k \V\(U v ,)L' 5; J k ), (B22) 



j„ = ^2j~H. (B23) 

Eq. l|BT3)l gives 

i) proton-proton or neutron-neutron matrix elements: 

((l^l^LS; J k \V$ kyTme \(lJ v ')L'S-, J k ) 
= Wo + c + (2S(S + 1) - 3)(6o + d )} ((l^)LS; J k \6(r - r')\(lJ u ,)L'S; J k ) 

+ { 0l + c x + (25(5 + 1) - 3)(&! + d x )} ((l^)LS; J k \tf 2 5(r - r') + S(r - r')k 2 \{l v l v ,)L' S; J k ) 
+ {a 2 +c 2 + (25(5 + 1) - 3) (6a + d 2 )} {(l^)LS; J k \tf ■ S(r - r')k\(U v >)L'S; J k ) 
+ {a 3 + c 3 + (25(5 + 1) - 3) (6 3 + d 3 )} {(l^)LS; J k \p^ Q (r)S(r - r')\(U„,)L'S; J k ) 

+2 e3 ((y M 0^;^bio(rK- 1 (r)«(r-r')|^^)i / 5;J fc ) x | ( ~| } ' 

+h((l fl l„')LS; J k \ P 2 10 (r)p^ 2 (r)S(r - r')\(l u l„,)L'S; J k ) 

+(a 4 + cJHy^LS; J k \{<r + a 1 ) ■ fct x 8{r - r')k\{l v l v ,)L' S; J k ), (B24) 

ii) proton-neutron matrix elements: 

((Ipln^LS; J k \V^ yrmc \(lJ u >)L' S; J k ) 
= {a - co + (25(5 + 1) - 3)(6o - doMl^LS; J k \S(r - r')\(lJ v ,)L'S; J k ) 

+ {at - c x + (25(5 + 1) - 3)(6 a - d x )} ((l^)LS; J k \k^5(r - r') + S(r - r')k 2 \{l v l v -)L' 'S; J k ) 
+ {a 2 -C2 + (25(5 + 1) - 3)(k - d 2 )} ((l^')LS; J k \tf ■ S(r - r')k\(lJ„,)L' S; J k ) 
+ {a 3 - c 3 + (25(5 + 1) - 3)(6 3 - d 3 )} ((l^)LS; J k \p^ (r)S(r - r')\{l v l v ,)L' S; J k ) 
+h{(y»>)LS; J k \pUr)p%o 2 {r)8{r - r')\{lJ„,)L' S; J k ) 

+(a 4 - c^dl^LS; J k \{<r + a') • fct x 5(r - r')k\(l v l v ,)L'S; J k ). (B25) 

We use the canonical (and real) radial wave functions i? M (r), the angular wave functions Yj^^fi), and the spin wave 
functions to write the nontrivial matrix elements included in Eqs. 1B24|) and i|B25[l as 

{(y^)LS; J k \S(r - r'^lJ^L'S; J k ) 

= J drr 2 R^r)R^(r)R u (r)R u ,(r)6 LL , J dQ JV, M (0)] * q \Y lv (Q)Y^, (fi)] l q , 

(B26) 
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dCl 



^(n)y, (n) 



LO 











(B27) 



((lyly.)LS; J k \S(r - r')k 2 \(l v l u ,)L' S- J k ) 



?LL> 



(in 



LO 



d 2 2d l v {l v + 1) 
dr 2 r dr r 2 



R v {r) 



R„>(r) 



-Ru{r) 



2d l v ,{l v , + l) 



dr 2 r dr 



R V ' (r) 



£ £ \^{l„ + M u ){l v ,+M v ,)l l »' 

AZ„=0,1 A2„/=0,l 



lu 



x / drr 2 R^r)R^(r)\(l v + l-Al v y-^ + (-)^ dRv{r) 



M r ) 



dr 



x^(^ + l-A^)^ + H Ai "'^^ 
r cir 



X<5 



LL' 



Y llt (si)Yt (a) [Y L ^ 1+2Alt/l (n)Y lu _ 1+2A an)] L0 , 



(B28) 



/o 

((lyly)LS; J fc |fct • J(r - r')fc|(l„WS; J fc ) = - WOmWIK*^ ■ <*(r - r'JfcJoIKU^)^, (B29) 



((lyly>)LS; J k \i(<r + <t') • fc* x <5(r - r'JfcKU^i'S; J fc ) 



(-) 



i+L'+./ fc4V 3| 1 ^' Jk \5 sl ((l l ,l fl ,)L\\(k^5(r-r')k) 1 \\(U ul )L'). (B30) 



The square brackets around products of spherical harmonics and the parentheses surrounding products of operators 
indicate angular-momentum coupling. 

To evaluate Eqs. (|B29I) and (|B30|) . one can use 



{{UM\{k^ ■8{r-r')k) I \\{l v l u ')L l ) 



\ E E + 

AZ U =0,1 Ai„=0,l i„, ,/ 



Mr) , , sAiJR^r) 



+ H 



x (Ul- AZ„) 



[ ' dr 



Rw{r) 



x + A/^V'i/ + A/^2/^ + 4AZ M - 1(> + 4A/„ - \l v , 



I L' L I I /„„/ 
V 1 1 



L 



Oflft' 

1^ l„ + 2Al tl -l Ip, 



V 1 

l v l v + 2/\l u - 1 l v , 





l v +l„,+L' 



(-) 





n and <-> v ], 

where for reduced matrix elements we have used the convention 

1 



v v'\- (-)'"+V+ i [ M «_> M ' 



(LL z |O im |L'L^) = T (L'L^m|LL z )(L||0 ; ||L'), 



(B31) 



(B32) 
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and made the abbreviation 

[A^ vv .] - (-) l »+ l *'+ L '[v «-> v'\ - «-> fjf] + [fx^fj,' and v v'\ 

— A , , — ( — \lv+K'+ L ' A , , — ( — \ l v-+ l n'+ L A , , J- A , , (TK^'W 

Eq. (|B26|) . modified to include additional factors in the radial integral, can also be used (together with the subsequent 
equations) to evaluate the matrix elements of the terms involving Poo{ r ) in Vgj^, rmc , the Coulomb-exchange interaction, 
and the contributions of the pairing functional to the effective ph, pp, and 3plh interactions. The Coulomb-direct 
term can be evaluated in a similar but slightly more complicated way, via a multipole expansion. 

In the main part of this paper we used the Skyrme functional SkM* , which is usually parameterized as in interaction 
in terms of coefficients to, ti,t2 j t^,x^ j x\,X2,x 3 , and Wq. The relations between these coefficients and those used here, 
if no terms are neglected, are j__3, Li} 

Co = |*o + Poo, | ^ Cf = -3*0(5 +%o)- i3*3(f + #3) P001 

Q> = ^4*0(5 _ ^o) - 24*3(2 _ x s) PoO' = ~8*0 ~ 48* 3 ^00' 

C ° = ll*i + 1*2(3 + *a), c l = -1*1(5 + *i) + 5*2(5 + 

Co= -5*1(5-^1) + 1*2(5 +X2), ^ = -1^*1+15*2, (B34) 

^0 S = 32*1(2 _ Xl ) + 32*2(2 + 2: 2), C x 8 — gjtl + gjt 2 , 

c w = ~|^o, cY J = -\w Q . 

In the HF fits that originally determined the SkM* parameters, the effects of Cj (the "J 2 terms") were neglected 
because of technical difficulties. These terms have often been included in subsequent RPA calculations. To maintain 
self consistency here, we have set them to zero both in the HFB calculation and in the QRPA. 
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